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Abstract 


A first principles dynamic model of the physical, chemical, and electrochemical processes at work in a proton exchange membrane fuel cell 
has been developed. The model solves the dynamic equations that govern the physics, chemistry and electrochemistry for time scales greater than 
about 10 ms. The dynamic equations are solved for a typical but simplified quasi-three dimensional geometric representation of a single cell repeat 
unit of a fuel cell stack. The current approach captures spatial and temporal variations in the important physics of heat transfer and water transport 
in a manner that is simple enough to make the model amenable to PEMFC system simulations and controls development. Comparisons of model 
results to experimental data indicate that the model can well predict steady state voltage—current relationships as well as the oxygen, water, and 
nitrogen spatial distribution within the fuel cell. In addition, the model gives dynamic insight into the distribution of current, water flux, species 
mole fractions, and temperatures within the fuel cell. Finally, a control system test is demonstrated using the simplified dynamic model. 


© 2006 Elsevier B.V. All rights reserved. 
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1. Background and introduction 


Due to quiet operation, high efficiency and low pollutant 
emissions performance characteristics, proton exchange mem- 
brane fuel cells (PEMFC) are receiving increased attention 
for applications in the automobile, portable power and backup 
power industries. PEMFC utilize air and hydrogen to generate 
electrical power that can be used to meet a wide variety of elec- 
trical load demands. Use of PEMFC in automotive applications 
especially, and in many other applications, will require a sig- 
nificant capability to follow dynamic load changes. In addition, 
to ensure safe and efficient operation and to develop control 
strategies for PEMFC systems it is important to understand both 
the steady state and dynamic performance characteristics of the 
PEMFC. 

Modeling has been used extensively to garner insight into 
PEMFC performance. Most of the PEMFC modeling efforts in 
the literature address steady state performance through multi- 
dimensional simulation. More recently, reduced order (bulk 
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or one-dimensional) transient fuel cell modeling and multi- 
dimensional fuel cell component transient modeling efforts have 
been presented in the literature. The current dynamic model is 
novel in that it attempts to resolve some of the geometric fea- 
tures of a PEMFC in a simplified manner sufficient to make 
it amenable to use in PEMFC system simulation and controls 
development. The geometric simplification employed in the cur- 
rent work is a 2+ 1D reduction (quasi-3D) with one dimension 
(1D) resolved through the membrane electrode assembly and 
two dimensions (2D) resolved in the fuel cell repeat unit cross- 
section accounting for gas and thermal transport in lumped 
serpentine channels. 

An excellent review of steady state and multi-dimensional 
PEMFC modeling has been developed presented by Weber and 
Newman [1]. Weber and Newman characterize these models by 
affiliation as being derived from the significant early contribu- 
tions of Springer et al. [2], Bernardi and Verbrugge [3,4], and 
those using computational fluid dynamics or other approaches. 
Since publication of the Weber and Newman paper, steady 
state PEMFC modeling has continued to be an active area of 
research as characterized by recent contributions from many 
groups including Lum and McGuirk [5], and Shimpalee et al. 
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Nomenclature 


BQ MR ESN 


Omo 


water activity, Effective catalyst area per unit vol- 
ume [m7 !] 

activation polarization coefficient 

surface area [m?] 

solid specific heat capacity [kJ kg7! K~!], capac- 
itance of the double-layer [F m2] 

species molar concentration (H2, H20, N2, O2) 
[kmol] m~3] 

constant volume gas specific heat capacity 
[kJ kg"! K7!] 

constant pressure gas specific heat capacity 
[kJ kg"! K7!] 

diffusion coefficient [m? s~!] 

hydraulic diameter [m] 

membrane dry equivalent weight [kg kmol~!] 
Faraday’s constant [96,485 C mol~!] 

mass transport coefficient [m s7!] 

change of Gibbs free energy [kJ mol 1] 
enthalpy [kJ kmol~!], convective heat transfer 
coefficient [kW m7? K7!] 

enthalpy of formation [kJ mol~!] 

exchange current density [A m~7] 

electrical current [A] 

fluid conduction heat transfer coefficient 
[kW m7! K7!] 

solid conduction heat transfer coefficient 
[kW m7! K7!] 

length [m] 

number of participating electrons in the reaction 
osmotic drag coefficient 

molar capacity, or total number of moles [kmol] 
species molar capacity [kmol] 

molar flow rate [kmol s~!] 

Nusselt number 

pressure [kPa] 

heat transfer [kW] 

universal gas constant (8.3145 kJ kmol~! K7!) 
Fuel cell external/load resistance [Q] 

species reaction rate (H2, H20, No, Or) 
[kmols~!] 

Sherwood number 

time [s], thickness [m] 

temperature (K] 

voltage [V] 

volume [m°] 

spatial distance [m] 

species mole fraction (H2, H20, N2, O2) 


reek letters 


species diffusion flux between GDL and bulk 
gases (H2, H20, N2, O2) [kmol s~!] 
water osmotic flux through the MEA [kmol s7!] 


Wmo water diffusion flux between GDL and MEA 
[kmol s7!] 

E porosity of GDL 

ô thickness [m] 

K ionic conductivity [S m!] 

À membrane hydration 

p density of solid [kg m~7] 

o electronic conductivity [S m!] 


Subcripts 

act activation polarization 
cl catalyst layer 

gdl gas diffusion layer 

in into control volume 

m membrane dry basis 
mea membrane electrolyte assembly 
mem membrane 

o standard condition 
ohm ohmic polarization 
out out of control volume 
sat water saturation 


[6], Pasaogullari and Wang [7], Guvelioglu and Stenger [8], 
Gorgun et al. [9], Nguyen et al. [10], and Freunberger et al. 
[11], among others. In addition, Burt et al. [12] and Campanari 
and Iora [13] have developed similar multi-dimensional mod- 
els for solid oxide fuel cells. While many additional groups are 
making significant contributions to steady state fuel cell mod- 
els, discussion of these contributions is limited herein since the 
current paper addresses transient PEMFC modeling. 

Amphlett et al. [14,15], Wohr et al. [16], and van Bus- 
sel et al. [17] were some of the first researchers to develop 
and use dynamic PEMFC models in the late 1990s. These 
dynamic models were generally bulk models that accounted 
for the dynamic physical, chemical and electrochemical pro- 
cesses in each single cell component (e.g., gas diffusion layer, 
membrane). These transient models were able to predict fuel 
cell voltage and heat losses as a function of time due to vari- 
ous changes imposed on the PEMFC. Subsequent reduced order 
(bulk or one-dimensional) dynamic PEMFC models have been 
developed by many groups. The bulk dynamic model of Yuyao 
and Choe [18] elucidated the mechanisms of PEMFC dehy- 
dration. Yerramalla et al. [19] developed one of the first bulk 
dynamic PEMFC models for control system development in 
Simulink®. Ceraolo et al. [20] and Pathapati et al. [21] each 
developed dynamic bulk PEMFC models accounting for the 
primary processes that govern membrane electrode assembly 
behavior in Simulink®. Pukrushpan et al. [22] were the first to 
demonstrate control system development based on the dynamic 
understanding of water transport in a bulk PEMFC model. Xue et 
al. [23] developed a dynamic PEMFC model that could predict 
the effects of temperature, gas flow, and capacitance on sys- 
tem transient behavior. Lemes et al. [24] used a bulk dynamic 
PEMFC model in a hardware-in-the-loop testing effort. In addi- 
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tion, several groups have developed bulk dynamic models for 
the related direct methanol fuel cell, including Sundmacher et 
al. [25], Kulikovsky [26] and Xu et al. [27]. While each of these 
models makes a significant contribution, none resolve the planar 
distribution of species, temperature, membrane water content, 
etc. within the PEMFC and some use simplified approaches for 
modeling the dynamics of water transport. 

More recently, multi-dimensional dynamic modeling of 
PEM fuel cells has been explored. The focus of these multi- 
dimensional dynamic models has been on resolving internal 
cell transient behavior. For example, Um and Wang [28] have 
investigated the interaction between mass transport and electro- 
chemical kinetics at the cell level. The more recent contribution 
of Um et al. [29] presents a transient, multi-dimensional model 
that accounts simultaneously for electrochemical kinetics, cur- 
rent distribution, hydrodynamics, and multi-component trans- 
port. Chen et al. [30] investigated the effect of water transport 
dynamics on cell transient capability, and Yan et al. [31] have 
investigated the effects of flow distribution and gas diffusion 
layer (GDL) morphology on the transient performance during 
startup operation. Berning and Djilali [32] presented a compu- 
tational fluid dynamics multiphase model of a PEM fuel cell 
that accounted for three-dimensional (3D) transport processes 
with phase change and heat transfer in the gas diffusion layers, 
the gas flow channels, and the cooling channels. Yu et al. [33] 
developed a dynamic water and thermal management model for 
a Ballard PEMFC stack. Shimpalee et al. [34,35] present a three- 
dimensional numerical simulation of the transient response of a 
PEM fuel cell (PEMFC) subjected to a variable load. Transient 
responses of the cell are predicted based on local distributions 
of the current density and gas mole fractions. The predictions 
show transients in the current density that overshoot the final 
state when the cell voltage is abruptly changed from 0.7 to 0.5 V 
for fixed initial flow rates in excess of stoichiometric to nor- 
mal [34]. Similar transient predictions are made for normal to 
minimal flow conditions in [35]. Wang and Wang [36] also devel- 
oped a three-dimensional transient PEMFC model, which when 
applied to a single channel PEMFC showed that the time for 
fuel cells to reach steady state was on the order of 10 s due to the 
effect of water accumulation in the membrane. Wang and Wang 
observed overshoot and undershoot of the current densities dur- 
ing step changes of some operating conditions. 

These multi-dimensional transient investigations have typi- 
cally utilized finite difference methods with detailed resolution 
of geometrical features in multiple cells and are quite compu- 
tational intensive. In addition, many are developed using com- 
mercially available computational fluid dynamics packages that 
make them not amenable or suitable for system level transient 
study and controls development. 

The objective of the current work is to develop a model that 
can capture both the dynamic response characteristics and some 
effects of cell geometry in a PEMFC model that is simplified 
enough to be useful in system level simulation and control sys- 
tem development. The present model is a balance between the 
high geometric resolution of prior steady state and dynamic cell 
and stack models and the low to zero geometric resolution of pre- 
vious reduced order dynamic models. The approach simplifies 


the geometry by only considering a quasi-three dimensional rep- 
resentation of a single fuel cell repeat unit of the fuel cell stack. 
In addition, only the physics, chemistry and electrochemistry 
extant in the fuel cell that contribute to dynamic performance 
of a PEM fuel cell at time scales greater than about 10 ms are 
resolved. All other processes are assumed to be in equilibrium 
or quasi-steady state. Finally, the model is developed in the 
Simulink® framework to allow development and testing of con- 
trol strategies and techniques. 

At the timescales of interest to the current effort (i.e., greater 
than 10 ms) the dynamic voltage response of a PEMFC depends 
primarily upon: (1) the amount of current drawn from the fuel 
cell, (2) localized species mole fractions, (3) local temperatures, 
(4) pressures at local triple phase boundaries, and (5) water con- 
tent of the membrane. Therefore, in order to garner detailed 
insights into dynamic fuel cell performance characteristics, it 
is necessary to understand and resolve the local behavior of 
species mole fractions, temperatures, pressures, and water trans- 
port within the fuel cell during transient operation. In addition, 
it is desirable to develop such a detailed dynamic model in a 
framework that allows the development and testing of control 
strategies for safety and enhanced performance. 

By discretizing the PEMFC in two dimensions and resolving 
heat and mass transfer and electrochemical processes in the third 
dimension (quasi-three dimensional) it is possible to predict the 
dynamic distribution of current, temperature, species mole frac- 
tions, and membrane water content in the fuel cell dynamically. 
The local resolution provides insight into the local performance 
of the fuel cell, local thermal stresses, local hydrogen and oxygen 
starvation, as well as potential flooding conditions and loca- 
tions during transient operation. In addition, the local resolution 
enables more accurate prediction of overall fuel cell dynamic 
performance compared to bulk dynamic models, which is bene- 
ficial for control and system design. To evaluate the model and 
verify model accuracy, experimental fuel cell data were acquired 
and compared to results from the current model. Once verified, 
the modeling technique as described herein can be applied to 
other PEMFC geometries and used to evaluate system config- 
uration and control strategies. To demonstrate the utility of the 
model for control design, a robust power and current-based fuel 
control is developed and demonstrated using the model. 


2. Model development 


The current model was developed to dynamically resolve 
local states and operating features of the fuel cell, but yet be 
simple enough for simulation on a personal computer (~2 GHZ) 
and incorporation into systems level models for control system 
development and testing. Since the model is to be utilized for 
controls development it was developed in Matlab-Simulink®. 
Simulink® is a preferred framework for controls development, 
but imposes a major limitation to the resolution that can be cap- 
tured in the model. As a result, neither all geometric features 
of the particular PEM fuel cell stack nor all of the physical 
and chemical processes occurring in the fuel cell stack can 
be captured without the model becoming too computationally 
intensive. Thus, many simplifying assumptions and geometric 
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Table 1 
Summary of order of magnitude time analysis 


Process Governing equation Numerical evaluation Time (s) 
1 1 1 1 
Charging or discharging of the d2,aC (G + 2) $ (1 x 1075 m” x 10° m7! 0.2Fm7? ( OS m H 30005m ) 2x 10-7 
electrochemical double layer eoe i oe 
Electrochemical reaction rate NA NA 1 x10773? 
5 2.5 x 10-4 my 
Species diffusion et seed 6 x 1073 
Ds 107° m? s~ 
ôm AA/EW 2200kg m~? x 85 x 1076 14/1000 kg kmol”! 
Water accumulation in the membrane Pmôm AXJEW li a ee i eons 34 
1/2F 15 A/2 x 96485 C mol™ 
PmdmCp/AT 2200kg m~? x 85 x 1076m x 2.2kJkg™! K7! x 10K 
Heat accumulation in membrane — 70 
ViostZ/A 0.4V x 15 A/0.1014m?/1000 J kJ- ! 
5eCp AT 2210kgm-3 x 5 x 1077 0.5kJkg~! K7! x 10K 
Heat transfer in electrodes pail haa ee a = 2 934 
Viostl/A 0.4V x 15A/0.1014 m2/1000JkJ~! 


à See Ref. [36]. 
b See Ref. [37]. 


simplifications have been made to enable solution of the dynamic 
equations that govern the PEMFC. 

Only processes that affect the larger timescales (>10 ms) 
associated with PEMFC performance are resolved dynamically. 
An order of magnitude analysis was accomplished in the cur- 
rent effort to determine whether processes must be dynami- 
cally resolved or can be assumed to be in equilibrium for the 
timescales of interest (>10 ms). The results of this order of mag- 
nitude analysis are presented in Table 1. Note that the complete 
kinetic analysis of electrochemical reactions would require the 
determination of the rate constants for each elementary reaction 
step. Hence a governing equation for electrochemical reaction 
rate is not presented in Table 1. However, the electrochemical 
time scale is on the order of 10~? s as reported by [36]. From the 
order of magnitude analysis, it has been determined that dynam- 
ics of electrochemical reactions, and charging or discharging of 
the charge double layer can be ignored since only time scales 
greater than 10 ms are of interest. On the other hand, species 
transport and diffusion, water accumulation in the membrane 
and heat transfer effects should be captured. Heat effects should 
be captured when investigating transients on the order of tens of 
seconds since heat generation in the membrane can have time 
scales in this range. 

The current model was developed in the Simulink® frame- 
work using the numerical method of lines as the basic solution 
technique. Each of the components that comprise the repeat unit 
of a PEMFC is discretized in space using control volumes. In this 
fashion the mechanisms of system energy and species conserva- 
tion can be described in terms of ordinary differential equations 
that can be solved in Simulink to determine local temperatures 
and species mole fractions in each control volume. From local 
temperatures and species mole fractions of each control volume, 
the reaction rate, molar capacity, heat transfer, and thermody- 
namic properties can be evaluated throughout the quasi-three 
dimensional cell. In addition, by assuming quasi-steady electro- 
chemistry (the kinetics of which only affects dynamic perfor- 
mance at timescales less than about 10 ms), the fuel cell voltage 
and current generation can be determined. Of course the losses 
associated with the electrochemical kinetics are accounted for 


through a bulk activation polarization term. In this fashion the 
model predicts local temperatures, species mole fractions, and 
flows as well as the extent of heat transfer and local electrical 
work generated throughout the cell. 


3. Assumptions 


The following is a list of the major assumptions made: 


1. Control volumes are characterized by a single lumped tem- 
perature, pressure, and set of species mole fractions condi- 
tion. 

2. All gases are ideal gases [11]. 

3. Auniform gas pressure is assumed. The pressure drop along 
the gas flow channels and GDL are neglected [11]. 

4. One-dimensional fully developed laminar flow occurs along 
the stream-wise direction [11]. 

5. Parallel diffusive fluxes in the gas diffusion layer and 
membrane are ignored. Convective transport inside the 
flow channel by far dominates parallel diffusive fluxes 
(Pe = 10!) [11]. 

6. The solid gas diffusion layer (GDL) and membrane elec- 
trode assembly (MEA) have a lumped temperature. (The 
respective Biot number was found to be much less than 
0.1.) 

7. Each cell in the stack is assumed to operate identically, so 
that a single cell simulation is taken as representative and 
used to calculate full stack performance [11,12]. 

8. All electrodes are good conductors for which an equipoten- 
tial electrode surface is assumed [13]. 

9. Quasi-steady electrochemistry is assumed, since the elec- 
trochemistry is rapid (occurring at time scales on the order 
of 107° s) [37]. 

10. A single activation polarization equation is used to capture 
the effects of all physical and chemical processes that polar- 
ize the charge transfer process. 

11. All reactants generate their ideal number of electrons, and 
no fuel or oxidant crosses the electrolyte [11]. 
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Flow Inlet 


Flow Outlet 


Anode Solid Plate 


Anode Bulk Gas 


Anode GDL 

Cathode GDL Flow Plane (x,y) 
Cathode Bulk Gas 
Cathode Solid Plate 


Flow Perpendicular Direction (z) 


Fig. 1. Illustration of the stream wise as well as the cross-sectional discretization 
of a PEMFC. 


4. Discretization 


To solve for local states of the fuel cell, the fuel cell is dis- 
cretized into control volumes quasi-three dimensionally. The 
fuel cell is discretized perpendicular to the flow through one 
repeat unit of the fuel cell stack as well as two dimensionally in 
the flow plane to capture the superficial area of the cell and the 
serpentine flow pattern as well. Flow perpendicular discretiza- 
tion is accomplished in a fashion similar to Yuyao and Choe [18] 
and Freunberger et al. [11]. The primary components of the fuel 
cell (bipolar plate, anode flow, anode GDL, MEA, cathode GDL, 
cathode flow, cathode plate and coolant flow) are discretized into 
perpendicular control volumes (Fig. 1) for each of the nodes used 
to discretize the fuel cell in the flow plane. In this fashion the 
fuel cell is discretized in the vertical direction into eight control 
volumes, using five types of control volumes: (1) solid plate, (2) 
bulk gas, (3) GDL, (4) MEA, and (5) coolant. 

Each of the eight control volumes in the flow perpendicular 
direction, were chosen because the operating voltage and local 
current production of a fuel cell is governed by the localized 
species mole fractions, temperature, and pressure at the triple 
phase boundary and the extent of membrane hydration, which 
can only be resolved with understanding of local heat trans- 
fer, electrochemistry, water transport, etc. in this perpendicular 
direction. In order to resolve the local species mole fractions at 
the triple phase boundary, the anode and cathode gas and GDL 
control volumes were desired. The gas control volume repre- 
sents the bulk flow of the gases. The gas from the bulk flow 
then diffuses through the GDL to the triple phase boundary, 
between the GDL and the MEA. Diffusion fluxes and osmotic 
water transport to the GDL control volume is determined for the 
boundary between the electrolyte and GDL control volume. In 
this fashion the species calculated to diffuse through the GDL 
control volume represent the species mole fractions at the triple 
phase boundary of the fuel cell. 

Itis possible to model a PEMFC without resolving the GDL or 
the triple phase boundary by assuming rapid diffusion between 
the triple phase boundary and the bulk flow. In this case, the 
concentration gradient between the bulk flow and the triple phase 


Table 2 
Summary of perpendicular discretization 


Perpendicular Species mass Energy conservation 
discretization conservation 
Solid plate Temperature Eq. (1) 
Bulk gas Molar flow rate Eq. Temperature Eq. (2) 
(7) species mole 
fraction Eq. (6) 
GDL Species mole number Temperature Eq. (5) 
Eq. (8) 
MEA Water content Eqs 
(11) and (12) 
Coolant Temperature Eq. (2) 


boundary as it impacts cell voltage can be somewhat accounted 
for using a bulk concentration polarization term [38]. Use of 
this assumption would make it possible to simplify the model 
by removing the GDL control volume. However, removing the 
GDL control volume makes it quite difficult to accurately model 
water transport and membrane hydration, which significantly 
affects conductivity and overall cell performance. Research has 
indicated that there can be significant water gradients between 
the bulk gas phase, triple phase boundaries and the electrolyte 
due to water formation at the cathode triple phase boundary, 
and osmotic drag [7]. As a result, to more accurately model the 
extent of membrane hydration it is beneficial to discretize the 
fuel cell utilizing both GDL and MEA control volumes. This 
makes it possible to predict water movement in the membrane, 
accounting for water formation, osmotic drag, diffusion, and 
water storage within the fuel cell. 

The solid plate, bulk gases, coolant and a lumped GDL and 
MEA control volume are needed to resolve the perpendicular 
temperature profile of the fuel cell (see Fig. 3). While a separate 
GDL and MEA control volume is essential to resolve the water 
transport and concentration gradients within the fuel cell, the 
GDL and MEA can be lumped together into one perpendicular 
control volume to resolve the perpendicular temperature profile. 
This is justified because the applicable Biot number of the GDL 
and MEA is much less than 0.1. The perpendicular discretization 
is summarized in Table 2. 

One set of eight control volumes in the perpendicular direc- 
tion represents a single node of the cross-sectional area of the 
fuel cell. The fuel cell cross-sectional area can then be resolved 
by solution of local performance characteristics in each node 
as it represents a discrete part of the fuel cell serpentine flow 
path. Individual nodes are therefore connected to one another 
in the primary flow direction in a two dimensional extension of 
the approach introduced by Yi and Nguyen [39] and used also 
by Freunberger et al. [11]. Appropriate accounting for heat and 
mass transport between and amongst adjacent nodes in the flow 
direction and across flow channels is made. The flow parameters 
(flow rate, species mole fractions, and temperature) in the bulk 
gas and coolant control volumes are then passed from node to 
node in the stream-wise direction. Since the flow paths through 
the typical PEMFC are not straight (e.g., serpentine), the current 
discretization results in two-dimensional resolution of the flow 
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field and fuel cell performance, and three-dimensional resolu- 
tion of the temperature distribution within the fuel cell repeat 
unit. Fig. 1 demonstrates the PEMFC discretization in the two 
directions for a two by two node representation of a fuel cell that 
contains a total of 32 control volumes. 

The particular experimental fuel cell that we desire to simu- 
late contains seven small channels for oxidant and fuel transport, 
each of which deliver gases and remove products in a seven pass 
serpentine fashion as shown in Fig. 4. To simulate this particu- 
lar geometric configuration of the fuel cell, the seven individual 
small channels are lumped together, into single flow channels 
in each of the nodes used in the model. In addition, the ser- 
pentine path is captured by the use of five discrete nodes along 
the width and seven discrete nodes along the length of fuel cell 
cross-section (Fig. 4). 

Thus the current cell model contains a five by seven nodal 
representation of the flow field and fuel cell performance char- 
acteristics for a total of 35 nodes, and 280 control volumes. The 
discretized model geometry is presented diagrammatically in 
Fig. 4. Cross-channel flow effects can be important in PEMFC 
operation. In the current model, however, no cross-channel flow 
is allowed between the serpentine paths. Note that the current 
model flow paths already lump several actual flow channels of 
the experiment into one larger channel. In effect, this accounts 
for some of the cross-channel flow effects amongst the channels 
that would otherwise have comprised the actual flow paths of 
the experiment, while resolving the primary overall direction of 
flow. 

The model was discretized in a manner similar to Freunberger 
et al. [11] with eight control volumes in the flow perpendicular 
direction and many nodes along the flow direction. The cur- 
rent model could benefit from a more detailed discretization in 
the flow direction, except for the goal of assuring amenability 
to system simulation and controls development. Detailed CFD 
models for transient or steady state analysis typically contain on 
the order of 10,000 states, while Freunberger et al. discretized 
their model into approximately 2000 states. Using 2000 states 
in the current model for system control design is unrealistic as 
this will increase the computational burden on Simulink by a 
factor of 10 compared to the present model, resulting in conver- 
gence difficulty, and long simulation time. Hence a simplified 
discretization is used to retain practicality as a control and sys- 
tem design tool. 


5. Conservation equations 


Each control volume temperature, species mole fractions, 
molar flow rate, and water content is solved from conservation 
first principles. Energy conservation and heat transfer equations 
are used to solve for temperatures, species conservation and mass 
transport equations are used to solve for species mole fractions, 
molar flow rates and water content. The following describes the 
set of dynamic conservation equations that is solved in each con- 
trol volume. Subsequent sections describe the flux of species, 
reactions, and heat transfer between and amongst control vol- 
umes as needed to solve each of the conservation equations. 


5.1. Energy conservation 


Temperatures throughout the model are determined from 
solution of the dynamic energy conservation equation. The tem- 
perature of each solid plate control volume is found from solving 
the following ordinary differential equation (ODE): 


dT . 
¥C—=)> Q; 
P dt 2 in (1) 
where Qin represents heat transfer to the control volume. In a 
similar fashion each of the bulk gas and coolant control volume 
temperatures is determined by solution of the dynamic energy 

conservation equation: 


dT S : : ; 
NCV- = y Ninhin = 5 Nouthout F 5 Qin F 5 Oratent 
(2) 


where Ninhin represent the enthalpy flux into the control volume, 
Nout/out the enthalpy flux out of the control volume, Oin the heat 
transfer to the control volume, and Orntent represent the latent 
heat of liquid water evaporating into the bulk gas or vice versa. 
The temperature of the control volume is taken to be the exit tem- 
perature of the control volume, and the inlet temperature of the 
control volume is that of the upstream (from the “flow” perspec- 
tive) node exit temperature. Note that the “upstream” node can 
be different for each of the solutions for anode gas, cathode gas, 
and coolant streams since flows may traverse through the cell in 
different manners. This feature must be carefully accounted for 
in the model. The molar flow rate of each species is determined 
from species conservation. It is important to mention that the 
enthalpy flux in the gas control volumes is that of the bulk gas 
flow, as well as the gas that diffuses to-from the GDL control 
volume. The enthalpy of the fluid or gas of Eq. (2) is determined 
as: 


E T 
h = X Cp dT 3 
3 (Sax i ) A 


where X are the component mass fractions and Cp are the tem- 
perature dependent specific heats at constant pressure for each 
species. The number of moles in the control volume of Eq. (2) 
can be determined from the ideal gas law: 


P% 


N=— 
RT (4) 


The temperature of each lumped GDL and MEA control volume 
is found by combining the gas and solid control volume conser- 
vation equations. In addition, irreversibilities associated with the 
electrochemical reactions are modeled as heat generated in this 
control volume as follows: 
dT . . 

(> (pvc), + > (NC), Fr a Ninhin E NoutMout en 

: r | V-I 

G40... +AH-—-——— 
>, in £ latent nF 1000 (5) 
where AH is the enthalpy of formation of water, Z the fuel cell 
current, and Vis the fuel cell voltage. From each conservation of 
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energy ODE it is possible to directly determine the temperature 
of each control volume, which is assumed to be the same as the 
control volume exit temperature. 


5.2. Species conservation 


As energy conservation is used to determine the temperature 
of each control volume, species mole numbers at each gas con- 
trol volume exit is determined from species mass conservation. 
Specifically, the exit mole number of each bulk gas control vol- 
ume is found from the following species conservation equation: 


d(NX) . = z " 
ap = NinXin — NowXout + 2 (6) 


The exit molar flow rate is determined from the total species 
conservation equation: 


Now = Nin =- E ED D 


In Egs. (6) and (7), @ is the species diffusion flux from adja- 
cent GDL control volume. In order to solve for the exit mole 
fraction, and molar flow rate, the control volume mole fractions 
are assumed to be that of the control volume exit condition. We 
apply this often dubbed “perfectly stirred” assumption to all gas 
and liquid control volumes in the model. 


Local electrochemical reactions rates for each species ( R) are 
modeled in the GDL control volume. The species mole number 
of each GDL control volume is determined as follows: 


GT b+ mot Omo +R ® 


where Who is the water diffusion flux from the adjacent MEA 
control volume, and @y,0 is the amount of water osmotic flux 
through the MEA (between the two GDL control volumes) 
within each node. From species mole numbers, the species mole 
fraction and concentrations within GDL control volumes can 
readily be determined: 


N 
3 
N 
V 


X= 


(9) 


2 


C= (10) 


Species conservation is further used to calculate the amount of 
water in the MEA 


dNy,0 
Tr =L in a1) 


from which membrane hydration (à) can be determined by: 
-Nro EW 
Vn Pm (12) 


ny 


where Vm is the membrane dry volume, EW the membrane 
dry equivalent weight, and pm is the membrane dry density. 
This local membrane water content, which significantly impacts 


membrane ionic conductivity, is critically important to resolve 
to accurately determine fuel cell voltage model. 

From species mass conservation ordinary differential equa- 
tions it is possible to determine the species mole number in each 
of the bulk gas streams, each GDL, and in the MEA control 
volume of a single node. In addition, with the perfectly stirred 
assumption these ODEs define exit molar flow rates in bulk gas 
control volumes. The species mole number can then be used to 
determine mole fractions, concentrations, and membrane water 
content as appropriate. In general, the current approach and set 
of assumptions allows solution of dynamic equations that can 
account for all species in each major subcomponent of the fuel 
cell repeat unit. In the present model four species are consid- 
ered, namely: hydrogen, oxygen, nitrogen, and water since the 
modeled fuel cell operates on pure hydrogen and air. 


6. Heat transfer 


As described, energy conservation equations are used to 
determine temperatures throughout the fuel cell (Table 2). 
Specifically Eq. (1) is used to determine solid plate temperatures, 
Eq. (2) is used to determine bulk gas and coolant temperatures, 
and Eq. (5) is used to determine temperatures of MEA and GDL 
control volumes. In each of the three equations, the extent of heat 
transfer between adjacent control volumes needs to be defined. 

Heat is generated in the fuel cell due to irreversibilities in 
the electrochemical reactions. Recall that the impact of these 
irreversibilities on energy conservation is accounted for in the 
solution of the bulk GDL and MEA thermal control volume as 
the difference between the total energy available from the global 
electrochemical reaction, 


H2 + 502 — H20 (13) 


and the amount of energy exiting the fuel cell as electricity 
(AH(1/nF) — (VI/1000)). This concept can be visualized in 
Fig. 2. In the perpendicular direction (see Fig. 3) the heat gener- 
ated in the GDL and MEA layer, is transferred in part to the bulk 
gas through convection heat transfer and by diffusion of species 
from the GDL to the bulk gas. In addition, some heat from the 
electrolyte and GDL is transferred to the solid plate through 
conduction. The heat in the solid plate is then transferred to the 
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Fig. 2. V-I curve illustrating the amount of heat generated in the fuel cell. 
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Fig. 3. Cross-sectional heat transfer network. 


coolant channel, through convection. The electrode (solid plate) 
is in thermal contact with the GDL, coolant channel, as well as 
the bulk gas. As a result heat can transfer between the bulk gas 
and the solid plate electrode. The heat transfer modeled in the 
perpendicular direction is illustrated in Fig. 3. Note in addition, 
that heat is transferred between and amongst adjacent nodes as 
indicated in the two nodes of Fig. 3. 

The amount of heat generated at each location in the fuel cell 
can vary depending on the amount of current generated locally 
in the fuel cell. This variation, combined with the effects of 
flow oriented convective heat transfer and variations in gas and 
coolant flow rates and properties can lead to significant tem- 
perature gradients in the MEA and solid plates of the fuel cell. 
These gradients are somewhat “smoothed” by in-plane conduc- 
tion heat transfer. Thus, in addition to capturing heat transfer in 
the perpendicular direction as shown in Fig. 3, conduction heat 
transfer between adjacent solid plates and natural convection at 
the edge of each plate is captured in the flow plane as shown 
in Fig. 4. Conduction heat transfer is not modeled in the MEA, 
because this layer is very thin. 

Throughout the model convection heat transfer between solid 
and gas nodes are determined from Newton’ law of cooling as 
follows: 


Ò = Ah(T, — T1) (14) 


The convection coefficient (h) is determined from the Nusselt 
number (Nup) provided by Incropera and Dewitt [40] as listed 
in Table 3: 
_ Nupk¢ 
De 


h 


(15) 


where kg is the fluid conduction heat transfer coefficient, and Dy 
is the hydraulic diameter. 

Conduction heat transfer throughout the model is determined 
from Fourier’s law: 
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Fig. 4. Illustration of the fuel cell flow plane, flow discretization, and axial heat 


transfer network. 


The thermodynamic quantities used in the model are presented 
in Table 3. 


7. Species diffusion osmotic drag and reactions 


Species diffusion (È, Wmo), osmotic drag (Op,0), and reac- 


tion rates ( R) in the fuel cell, presented in the mass conservation 
equations are resolved in the model as presented in Fig. 5. The 


Table 3 
Model parameter values 
Description Value Units 
Geometry 
Cell width (x) 0.0507 m 
Cell height (y) 0.0507 m 
Depth of anode gas channel (z) 0.001 m 
Depth of cathode gas channel (z) 0.001 m 
Depth of cooling channel (z) 0.003 m 
Thickness of GDL (z) 0.001 m 
Thickness of electrolyte (z) 7x107 m 
Thickness of separator plates (z) 0.002 m 
GDL porosity 0.5 
Thermodynamic properties 
Separator place density 2210 kg m~? 
Separator plate specific heat capacity 0.5 kJ kg! K7! 
Electrolyte dry density 2200 kg m~? 
Electrolyte dry equivalent weight 1000 kgkmol~! 
Electrolyte solid specific heat capacity 2.179 kJ kg! K7! 
Heat transfer properties 
Separator plate conduction coefficient 0.22 kW m~! K7! 
Nusselt number of anode gas 2.96 
Nusselt number of cathode gas 2.96 
Nusselt number of coolant liquid 7.54 


Polarization constants (manipulated to match model to experiment) 
Exchange current density (io) 1.85 Am? 
Ohmic lost constant (b2) 350 
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Fig. 5. Cross-sectional species diffusion, species reactions, and osmotic drag. 


electrochemical reactions are modeled in the GDL control vol- 
umes, representing the region of triple phase boundaries in the 
current fuel cell model. The anode half reaction is: 


Hə > 2Ht 4+ 2e7 (17) 
and the cathode half reaction is: 
50) + 2H+ + 2e7 + H20 (18) 


with a global reaction as shown in Eq. (13). From Faraday’s law 
the reaction rate of both half reactions is directly proportional 
to the current as follows: 


Rmo = — Ru, = -2Ro, = —= (19) 

Species diffusion is captured in the perpendicular direction 
between gas, GDL, and MEA control volumes. Species transport 
from the gas channel to the GDL accounts for the convection 
driven by a concentration gradient and diffusion in the GDL. 
The mass transport coefficient (gm) at the gas channel and GDL 
interface is obtained based on the Reynolds analogy between 
heat and mass transfer: 


= SE Da 
Em = Du 


(20) 


where Sh is the Sherwood number, Din the diffusion coefficient, 
and Dy is the hydraulic diameter of the gas flow channel. The 
diffusion coefficients for species are functions of temperature 
and pressure and are modified via the Bruggeman correlation to 
account for the effects of porosity and tortuosity in the GDL as 
follows: 


>- a (TV IB 

Da = b(Z) (2) (21) 
o 

DË = e!5 Dm (22) 


where Do is the species diffusion coefficient at standard pressure 
and temperature, De the effective species diffusion coefficient, 
and € is the GDL porosity. The species diffusion flux between 


the GDL and bulk gasses is then as: 
1 
A—~ = 
(1/8m) + (tgai/ D) 


Rait = (23) 


@ = Rat(C2 — C1) (24) 


where Raif is the total diffusion resistance of each species and 
tga] is the thickness of the GDL. 


7.1. Water transport 


Since water content in the membrane strongly affects ionic 
conductivity, the current dynamic model captures the details of 
water behavior in the MEA. Two types of water molecule trans- 
port from anode or cathode GDL to the MEA are considered: (1) 
the electro-osmotic drag, and (2) diffusion due to a concentration 
gradient between control volumes.The rate of water molecule 
transport via osmotic drag from anode-—electrolyte interface to 
cathode—electrolyte interface is proportional to current density 
and the electro-osmotic drag coefficient as follows: 


1 
Omo = na — 25 
H0 = nig p (25) 


The osmotic drag coefficient (na), is calculated from the mem- 
brane water content (A), which depends on the water activity (a), 
as follows [41]: 


na = 0.002927 + 0.05 — 3.4 x 107!9 (26) 


for0 <a< 1 


forl<a<3 


0.043 + 17.8la — 39.85a2 + 36a?, 
144+ 1.4(a— 1), 


The water diffusion due to the concentration gradient between 
the two GDL and MEA control volumes is calculated by: 
C2— C1 
Wmo = DwA———— (27) 
tmea 
where Dw is the diffusion coefficient of water in the electrolyte 
and tinea 1S the thickness of the MEA. The water concentration in 
the membrane is calculated from the membrane water content: 


Cote (28) 


The diffusion coefficient of water in the electrolyte (Dw) is cal- 
culated from the empirical equation [41]: 


1 1 
Dy =D Wig) — = — 29 
w= Deap (= 7) (23) 
where 
1076, forà < 2 
3 [1+ 2(a — 2)] x 10-, for2 <1’ <3 
+) B= 1.670. —3)] x 1076, for3 <2 < 4.5 
1.25 x 1076, for4.5 <A 
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8. Electrochemical model 


By resolving species mole fractions, molar flow rates, and 
temperatures dynamically throughout the fuel cell and by 
assuming the fuel cell voltage is in quasi-equilibrium with the 
dynamic state of the fuel cell (Assumption 9), it is possible 
to determine the fuel cell operating voltage, and current gen- 
eration distribution throughout the fuel cell. It is important to 
note that because the fuel cell electrodes are good conduc- 
tors the voltage difference across any one cell is assumed to 
be constant for all parts of the cell (equipotential Assump- 
tion 8). As a result each nodal current must be determined 
such that all node voltages are equivalent. Before explaining 
the solution procedure, the voltage—current relationships are 
explained. 

Each nodal voltage is determined by subtracting locally cal- 
culated activation, and ohmic polarization from the Nernst volt- 
age: 


Vnode = VNernst = Vact a Vohm (30) 


The Nernst voltage at each node is determined based on the MEA 
temperature, and the local species mole fractions in the GDL 
control volumes, which are representative of concentrations at 


the triple phase boundary: 
1/2 
XH, X 
2 0 pi ) (31) 
gdl 


Xmo 


AG(T) RT 
VNernst = nF + nF In 


The activation polarization is modeled from the Tafel equation 
based on the local GDL control volume states: 


RT, I 
Vaat = ae =" In ( / ) (32) 
nF 


lo 


where aec is the activation polarization coefficient, and io is the 
current exchange density. The ohmic polarization is modeled as 
determined in [2] for Nafion 117 based on the electrolyte control 
volume temperature and hydration as follows: 


V, h — I tmea 
om by exp(b2((1/303) — (1/Tmea))’ 
bı = 0.005139) — 0.00326 (33) 


The voltage is determined at each node from Eqs. (24)-(27), 
but to satisfy the equipotential assumption, each nodal voltage 
must be equivalent. In addition, the fuel cell must satisfy Ohm’s 
law for the external circuit: 


Vnode = Veell = 5 Thode R 84 


That is, the sum of all the nodal currents multiplied by the 
external resistance must be equal to each nodal voltage, or equiv- 
alently the cell voltage. As a result, each nodal current is iterated 
until all the nodal voltages are equivalent and ohms law is sat- 
isfied. It is important to note that the amount of current at each 
node effects species mole fractions, and temperatures throughout 
the fuel cell, which in turn affect the voltage. Fig. 6 illustrates 
the algebraic loop that must be solved to resolve the fuel cell 
voltage and current distribution. While this solution strategy is 
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Fig. 6. Voltage solution procedure. 
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being executed for each time-step, the external resistance can 
be manipulated, using a feedback loop, to control the fuel cell 
current, power, or voltage. 


9. Experimental setup 


A unit cell of a PEM fuel cell with an active area of 
25cm? was used to validate the simulation results. The unit 
cell contained 7-channel serpentine flow channels for gas deliv- 
ery to both the cathode and anode of the cell with a counter- 
flow direction. This basic configuration is presented schemat- 
ically in Fig. 4. The channel depth is 1mm and width is 
0.7 mm, and the channel turns are 7mm wide near the cell 
edge. The MEA used in this study is based on Nafion 112. 
The detailed specifications of the experimental cell configura- 
tion are presented in Table 3. Fig. 7 shows a schematic dia- 
gram of experimental setup, which consists of a gas supply 
unit, a unit fuel cell, an electronic load device, various sen- 
sors, a personal-computer (PC) based data acquisition system, 
gas sampling equipment, and a gas chromatograph (HP 5890 
Series II). 

High-purity hydrogen at the anode and dry air at the cath- 
ode are used as reactant gases and are humidified by passing 
each gas stream through a bubble type external humidifier. The 
flow rates of hydrogen and air were controlled and measured by 
two mass flow controllers. The humidification temperatures of 
the reactant gases were controlled by adjusting the humidifier 
temperature. 

In order to prevent water from condensing in the sample line 
for species measurement, a line heater with controller was used 
to maintain the temperature of the gas sample above 100°C. 
Species sampling points on the cathode side of the fuel cell are 
shown in Fig. 8. The size of each sample hole is 0.7 mm, which 
is the same as the width of the flow passage. The sampling points 
are evenly located along the gas channel spaced at 1/8th of the 
total channel length. To minimize the effects of gas sampling 
on the cell performance, the gas sampling rate was less than 5% 
of the total flow rate in the cathode channel. A cartridge heater 
was used to maintain the cell temperature as a constant dur- 
ing the experiment, which was controlled by a PID temperature 
controller. 
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Fig. 7. Schematic diagram of test equipment used for gas chromatograph measurements. 


10. Steady state validation and analysis oxygen, nitrogen, and water vapor molar concentration along 
the length of the cathode gas channel. The polarization curve 

Comparisons between the simulation and the experiment comparison result is shown in Fig. 9. The model well simulates 
were made in terms of the cell polarization curve as well as the overall cell polarization. Note that this good comparison was 
achieved by selecting a single point (labeled “reference point”) 
for establishing the constants in the cell polarization Eqs. (34) 
and (35) and holding them constant for all other conditions sim- 
ulated. Relative humidity and temperature of both anode and 
cathode inlets are held at 96% and 70°C, respectively, and the 
unit cell is also maintained at a constant operating temperature of 
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Fig. 8. Schematic diagram of gas sampling port positions. of 70°C, and relative humidity of 96%. 


F. Mueller et al. / Journal of Power Sources 163 (2007) 814-829 825 


0.7 
| CHeeees 
- fi ttteeteeeesesenesehosesheeeo 
= 4 
20.54 
= 
oO 4 
g 
w 0.44 
2 4 
£03 
f ] yyyvyvyyYYYY 
eyyy 
0.24 we 
ext 
| WANES 440465 E 
o bbe EPPP] 


SF 
0.0 0.1 02 03 04 05 06 0.7 08 09 1.0 1.1 
Fractional Distance 


Fig. 10. Comparison of mole fraction between experiment and simulation. 


70°C. Air and hydrogen utilizations are 0.33 and 0.53, respec- 
tively. Instead of modeling the cartridge heater of the experiment, 
a coolant channel is modeled as part of the cell support assem- 
bly (bipolar plate). The coolant fluid temperature of the model 
is adjusted to maintain the same constant operating temperature 
as that of the experiment. 

Comparison of oxygen, nitrogen, and water vapor mole frac- 
tion between experiment and simulation at different sampling 
locations is shown in Fig. 10. The relative humidity and temper- 
ature of both anode and cathode is 50% and 70 °C, respectively, 
in this case. The operating current is 15 A. Oxygen mole frac- 
tion at the cathode inlet is 0.167 and monotonically decreases 
along the cathode flow channel due to electrochemical reac- 
tion in the cathode compartment. Along the flow channel some 
oxygen is consumed at the MEA and two H20 molecules are 
produced for oxygen molecule consumed. Thus, the total molar 
flow rate increases and since water vapor is not condensed in the 
channel due to low relative humidity the oxygen and nitrogen 
mole fraction along the cathode channel is decreases consider- 
ably. If there were no H20 transport between the GDL layers by 
electro-osmotic drag or by back diffusion, the theoretical value 
of oxygen mole fraction at the cathode exit would be 0.113. Note 
that this value is close to the simulation result at the cathode out- 
let. This indicates that net H20 transferred from the anode GDL 
to cathode GDL is very small compared to the amount of H2O 
generated in the cathode due to electrochemical reaction. As a 
result, the magnitude of H20 transport by osmotic drag is almost 
equivalent to that of back diffusion from cathode to anode in the 
current experiment. The simulation result of oxygen, nitrogen, 
and water vapor molar concentration along the cathode channel 
well predicts the experimental data. 

One of the main features of the current model is the abil- 
ity to predict the local characteristics of current, temperature, 
and species concentrations throughout a PEM fuel cell without 
a complicated computational fluid dynamics (CFD) simulation. 
Because of the simplified resolution, the model can be integrated 
into system level models. To demonstrate the model capability 
to resolve current density distribution, membrane water con- 
tent, and water flux providing insight into potential operating 
issues, a high current density operation (1 A cm7? at 0.4 V) was 
simulated. The fuel cell operated at 70°C, with the inlet cath- 
ode and anode gas humidified to 60% at 70°C. The current 
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Fig. 11. Distribution of (a) current density and (b) membrane water content at 
pressure of 1.1 bar, cell temperature of 70°C, and relative humidity of 60%. 


density and water content distribution are presented in Fig. 11. 
The local current density is fairly uniform but increases slightly 
along the anode gas stream channel. As Fig. 11b indicates, the 
membrane water content is also quite uniform, so that the MEA 
is adequately hydrated over the entire active area even though 
the relative humidity of both the anode and cathode inlet flows 
is low. Thus, the local ohmic loss in the membrane, which is 
strongly dependent upon the water content, is almost constant 
and the high current density observed near the anode outlet is 
due to high oxygen concentration (affecting the Nernst potential) 
near the cathode inlet. 

Local distribution of hydrogen is shown in Fig. 12. Hydrogen 
mole fraction along the length of the anode channel is almost 
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Fig. 12. Hydrogen mole fraction in anode channel. 
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Fig. 13. The water flux magnitude of electro-osmotic drag, diffusion driven by 
the concentration gradient, and net flux from the anode to MEA at pressure of 
1.1 bar, cell temperature of 70°C, and relative humidity of 60%. 


constant, even though the overall hydrogen utilization is 0.52. 
This is due an almost constant net water flux between the anode 
to the cathode. About 28% of the H20 that entered via the 
anode gas stream was transferred to the cathode side due to 
high electro-osmotic drag, even though back diffusion trans- 
ported some water back to the anode GDL. As a result, even 
though more than half of the hydrogen is consumed, H20 mole 
fraction in the anode gas stream increases only slightly from 
0.167 to 0.25. Fig. 13 shows the water flux magnitude of the 
electro-osmotic drag process, the diffusion process as driven by 
the concentration gradient, and the net flux from the anode to 
the MEA. As the water flux due to the electro-osmotic drag is 
proportional to the protonic flux, the water flux due to electro- 
osmotic drag gradually increases along the anode channel length. 
But the water flux due to back diffusion peaks at the anode inlet 
and decreases along the anode channel length. As a result, the 
net flux between anode GDL and cathode GDL changes sign 
at a distance of about one-fourth of the channel length, which 
shows the internal circulation of water that well hydrates the 
membrane. 

Water molar concentrations at various locations in the cell 
components (perpendicular direction) are shown in Fig. 14 for 
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Fig. 14. Water molar concentrations at various locations in the through cross- 
sectional direction pressure of 1.1 bar, cell temperature of 70°C, and relative 
humidity of 60%. 


several specific nodes (numbers 1, 18, and 35). Note that the 
water concentration in the membrane is the fictitious vapor con- 
centration that is in thermodynamic equilibrium with membrane 
water content. One-to-one correspondence between the fictitious 
vapor concentration and water concentration based on water 
content can be obtained from Eqs. (27) and (29). The water con- 
centration gradient is significant across the membrane because 
the mass diffusivity of water in the membrane is two orders 
of magnitude smaller than that in the other regions. Therefore, 
the majority of the water that is produced at the cathode mem- 
brane interface is transported through the cathode GDL into the 
cathode gas channel. At the current temperature of operation, 
the water vapor is already saturated in the cathode GDL by the 
middle of the channel length. Although water content in both 
the anode and cathode compartments increases, the water molar 
concentration along the cathode channel length (from node 35 to 
node 1) increases more significantly than it does along the anode 
channel length (from node 1 to node 35). On the anode side the 
water vapor is transferred from the GDL to the channel near the 
anode inlet and then the direction of water transfer is reversed 
due to the substantial hydrogen reduction along the length of the 
anode channel and back diffusion. But on the cathode side water 
vapor is always transferred from the GDL to the gas channel. 

Oxygen mole fraction is substantially decreased along the 
air flow direction from the cathode inlet due to cathodic elec- 
trochemical reactions that consume oxygen and produce H20. 
Thus, in this high current density case with good membrane 
hydration, oxygen mole fraction ends up being the major fac- 
tor that determines the local current density. Relative humidity 
in the cathode gas stream is above 100% near the cathode out- 
let, which indicates that water condensation would likely occur 
under these operating conditions causing flooding in the cath- 
ode channel. The current model indicates that the majority of the 
water produced on the cathode is transported through the cathode 
GDL into the cathode gas channel. Cathode flooding by liquid 
water would significantly hinder the access of oxygen to the cath- 
ode catalyst layer and reduces the number of electrochemically 
active sites, resulting in a substantial decrease in local current. 
The current model does not contain a detailed understanding 
of water condensation physics or droplet formation, but, it can 
indicate the likely location of channel flooding as near the point 
of water saturation. No mechanism for affecting the cell perfor- 
mance due to flooding was included in the current model. As a 
result, the predicted local current density near the cathode out- 
let (beyond the point where water saturation occurred) may be 
higher than observed current densities. 


11. Transient comparison and controls 


While the model can be used for steady state analysis, the 
main feature of the model is transient capabilities resolving a 
simplified geometry in such a way it can be integrated into sys- 
tem level models. Results showing system level modeling the 
quasi-three dimensional PEM model are not presented herein 
because the focus of the paper is to present the quasi-three 
dimensional model and validate it. The models voltage response 
to an instantaneous 10-15 A increase was compared to the single 
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Fig. 15. Comparison of the voltage response between experiment and simulation 
during an instantaneous 10-15 A current increase. 


cell experiment (Fig. 15). The fuel cell operated at 70°C with gas 
at a 50% relative humidity. To avoid hydrogen and oxygen star- 
vation the utilization of each was 0.4748 and 0.2891 respectively. 
Following the current increase the voltage dropped from 0.69 to 
0.58 V instantaneously. The instantaneous decrease in voltage 
is a result of increased ohmic and activation polarization. Fol- 
lowing the instantaneous decrease the voltage slowly increased 
to 0.62 V in about 75s. The increase in voltage is a result of 
increased membrane hydration which results in a decrease in 
the cells internal resistance. As the current increased, the rate 
of water formation increases at the cathode which results in an 
increase in membrane hydration. The simulated transient was 
on the same order of magnitude as expected from the simple 
transient analysis presented in Table 1. Fig. 16 shows the mem- 
brane hydration during the voltage undershoot and at steady state 
following the current increase, illustrating the increase in mem- 
brane hydration that causes the slow voltage transient. Overall 
the simulated response represents very well the experimental 
response. There is a very slight offset in the voltage steady 
state response and a slightly larger undershoot in the simulated 
transient, but overall the simulated responses well match the 
experiment. 

During load changes, transients in membrane hydration are 
impossible to control exactly. Therefore, robust control strate- 
gies will have to be implemented in fuel cell systems to reject 
disturbances such as small local variations in membrane hydra- 
tion. The model developed herein is a powerful tool that can be 
used for such controls development. To demonstrate this fea- 
ture of the current model a power controller supplemented by 
current-based fuel control is implemented using the cell model 
in Simulink®. The power controller is shown in Fig. 17. Current 
is manipulated based on power using a proportional and integral 
feedback together with feed forward on the power demand. Inte- 
gral feedback makes it possible to track power with zero steady 
state error. 

Controlling the fuel flow rate in proportion to current makes 
it possible to operate the fuel cell at a desired high utilization, 
increasing the efficiency of the fuel cell. Utilization is the ratio 
of hydrogen consumed in the fuel cell to the amount of hydrogen 


Fuel 
Inlet, 


5 


Fig. 16. Water content in the membrane at the undershoot state and steady state 
at 15A. 


into the fuel cell: 


U= H2 consumed 3 5) 
Ain 


From Eq. (19), it is known that the amount of hydrogen con- 
sumed is proportional to the current. Hence utilization can be 
represented as: 


i/nF cells 
U = —— (36) 
(N Xn )in 
current 
l power e power d current + U current 


Y power 


Fig. 17. Robust power controller. 
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Fig. 18. Controlled load increase from 7 to 10 W. 


Therefore, to operate the fuel cell at constant utilization, the flow 
rate can be controlled by: 


i 


ÑN = —— (37) 
Un FXp, 


Eq. (37) is a simple rearrangement of Eq. (36) for a fixed utiliza- 
tion. Current-based fuel control, as shown in Eq. (37), together 
with the robust power controller shown in Fig. 17 were imple- 
mented to control a power demand increase using the quasi-three 
dimensional PEM model. A controlled power increase from 7 to 
10 W is demonstrated as shown in Fig. 18. The set-point power 
demand was tracked almost exactly by the model. During the 
transient the utilization remained nearly constant, with only a 
slight variation in the fuel cell. 


12. Conclusions 


A quasi-three dimensional dynamic model of a PEM fuel cell 
has been developed in Matlab-Simulink®. The model is based 
upon the dynamic equations that govern the first principles of 
mass and energy conservation, mass and heat transfer, and elec- 
trochemical reaction with special attention to the details of water 
transport within the cell. In addition, the fuel cell model captures 
some of the geometric features of the fuel cell as these dynamic 
equations are solved in a simple characteristic geometric config- 
uration. The fuel cell is discretized into (1) 5—7 control volumes 
through the membrane electrode assembly, and (2) 35 nodes in 
the stream-wise direction. Five perpendicular control volumes 
are used to solve the dynamic species and mass conservation 
equations and six perpendicular control volumes are used to 
solve the dynamic energy balance and to capture the details of 
MEA/GDL behavior, such as water transport, which is critical to 
accurately determine polarization losses and performance. The 
dynamic conservation equations, primary heat transfer equations 
and equations of state are solved in each control volume of the 
node, and each of these nodes is interconnected with adjacent 
nodes (in the stream-wise direction) to simulate local cell per- 
formance characteristics throughout the cell cross-section. 


Experiments were conducted to measure the performance of 
a single cell PEMFC under laboratory operating conditions for 
comparison to simulated results. A comparison of simulated and 
observed polarization curves shows that the dynamic model well 
predicts overall cell performance. In addition, since the model 
is discretized it can predict distributions of local current density, 
membrane water concentration, and species concentrations at 
both the anode and the cathode. Comparison of the simulated 
and observed oxygen water and nitrogen concentrations in the 
cathode compartment shows that the model well predicts species 
concentrations in the PEM fuel cell. The ability of the model 
to predict dynamic variations of performance characteristics in 
time and space based on first principles can be useful in the 
design of a fuel cell stack and for optimizing and controlling the 
dynamic cell performance for a variety of operating conditions. 
Also this model is shown to be a useful tool for investigating 
the effects of inlet conditions, operating parameters, and overall 
load and temperature conditions on both overall and local per- 
formance characteristics. Finally, such a dynamic model may be 
useful for the development of control strategies for enhancing 
overall PEM system performance. 

The current model can capture the local current density, mem- 
brane water content, and species distributions as they depend 
upon operating conditions, which can be utilized to improve 
cell design and optimize performance for various operating 
conditions. This model can be used to analyze the character- 
istics of the individual fuel cell transients and provide the basis 
for designing control strategies for integrated fuel cell systems. 
The first principles approach provides fundamental insights into 
local performance characteristics and allows capturing of the 
important physics, chemistry, and electrochemistry that governs 
the dynamic performance of a PEM fuel cell. The simplified 
geometry approach provides insight into spatial variations and 
challenges while at the same time being simple enough to include 
in simulations for the development of control systems. 
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